/*

  Simple cross-product, outputs in correct triangular form.

  - Coalescing memory access in the Column reads
  - No memory coalescing in the Row reads
  - 50% coalescing in the write (float8s)
  - All threads read in all data

 */

#define LOAD_ONE_BY_ONE_FLOAT(s)		\
  {*(input##s##_p) = *array_p;			\
    CUBE_ADD_BYTES(4);}

#define ONE_BY_ONE_COMPUTE(s)				\
  {float colXreal = input[s][4*tx + 0];			\
  float colXimag = input[s][4*tx + 1];			\
  float colYreal = input[s][4*tx + 2];			\
  float colYimag = input[s][4*tx + 3];			\
  float rowXreal = input[s][4*ty + 0 + 4*TILE_WIDTH];	\
  float rowXimag = input[s][4*ty + 1 + 4*TILE_WIDTH];	\
  float rowYreal = input[s][4*ty + 2 + 4*TILE_WIDTH];	\
  float rowYimag = input[s][4*ty + 3 + 4*TILE_WIDTH];	\
  sumXXreal += rowXreal * colXreal;			\
  sumXXreal += rowXimag * colXimag;			\
  sumXXimag += rowXimag * colXreal;			\
  sumXXimag -= rowXreal * colXimag;			\
  sumXYreal += rowXreal * colYreal;			\
  sumXYreal += rowXimag * colYimag;			\
  sumXYimag += rowXimag * colYreal;			\
  sumXYimag -= rowXreal * colYimag;			\
  sumYXreal += rowYreal * colXreal;			\
  sumYXreal += rowYimag * colXimag;			\
  sumYXimag += rowYimag * colXreal;			\
  sumYXimag -= rowYreal * colXimag;			\
  sumYYreal += rowYreal * colYreal;			\
  sumYYreal += rowYimag * colYimag;			\
  sumYYimag += rowYimag * colYreal;			\
  sumYYimag -= rowYreal * colYimag;}


CUBE_KERNEL(shared1x1float, cuFloatComplex *array_d, cuFloatComplex *product_d, 
	    const int write)
{
  CUBE_START;

  //set frequency number from blockIdx.y
  int f = blockIdx.y;

  int Row, Col, blockX, blockY;
  CUBE_DEVICE_CALL(findPosition, Col, Row, blockX, blockY);

  //instantiate sum variables
  float sumXXreal = 0.0, sumXXimag = 0.0;
  float sumXYreal = 0.0, sumXYimag = 0.0;
  float sumYXreal = 0.0, sumYXimag = 0.0;
  float sumYYreal = 0.0, sumYYimag = 0.0;

  //get local thread ID
  int ty = threadIdx.y;
  int tx = threadIdx.x;
  int tid = ty*TILE_WIDTH + tx;

  //declare shared memory for input coalescing
  __shared__ float input[2][TILE_WIDTH*8];

  float *array_p = (float*)array_d + tid + f*NSTATION*NPOL*2;
  if (tid < 4*TILE_WIDTH) {
    array_p += blockX*TILE_WIDTH*NPOL*2;
  } else {
    array_p += blockY*TILE_WIDTH*NPOL*2 - 4*TILE_WIDTH;    
  }

  float *input0_p = input[0] + tid;
  float *input1_p = input[1] + tid;

  LOAD_ONE_BY_ONE_FLOAT(0);
  array_p += NFREQUENCY*NSTATION*NPOL*2;

  //#pragma unroll
  for(int t=0; t<NTIME-2; t+=2){
    __syncthreads();

    ONE_BY_ONE_COMPUTE(0);
    LOAD_ONE_BY_ONE_FLOAT(1);
    array_p += NFREQUENCY*NSTATION*NPOL*2;
    
    __syncthreads();

    ONE_BY_ONE_COMPUTE(1);
    LOAD_ONE_BY_ONE_FLOAT(0);
    array_p += NFREQUENCY*NSTATION*NPOL*2;
  }

  __syncthreads();

  ONE_BY_ONE_COMPUTE(0);
  LOAD_ONE_BY_ONE_FLOAT(1);
  array_p += NFREQUENCY*NSTATION*NPOL*2;

  __syncthreads();

  ONE_BY_ONE_COMPUTE(1);

  if (write)
    CUBE_DEVICE_CALL(write1x1, Col, Row, product_d,  
		     sumXXreal,  sumXXimag, sumXYreal,  sumXYimag,  
		     sumYXreal, sumYXimag,  sumYYreal,  sumYYimag);

  if (Col <= Row) {
    CUBE_ADD_BYTES(32);         //output bytes
    CUBE_ADD_FLOPS(NTIME*32);
  }

  CUBE_END;
}
